Data from 1960 to 2010 raining days, there are 4168 observations 13 days of which are rain exceeds 1000 (10cm)
Change of exceeds 1000 is 0.003119002 0.31% change
library(tidyverse)
library(readr)
library(drc)
library(reshape2)
#install.packages("moments")
library("moments") # for skewness of data
Data
read.csv("./raw/prcp.csv") %>%
mutate(date = as.Date(date)) %>%
mutate(year = format(date, "%Y")) %>%
filter(prcp != 0) %>%
dplyr::select(date, year, prcp) -> data
data %>%
filter(year <= 2020) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
Modling Distribution
Use data from 1960 - 2010 for modeling precipitation
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
summary(model)
Model fitted: Exponential decay with lower limit at 0 (2 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
d:(Intercept) 0.1407885 0.0012481 112.80 < 2.2e-16 ***
e:(Intercept) 1.3292414 0.0077504 171.51 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error:
0.001397153 (945 degrees of freedom)
plot(x = dens$x, y = fitted(model), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

mod.exd <- drm(y ~ x, fct = EXD.2(), data = dens)
plot(x = dens$x, y = fitted(mod.exd), type = "l", col = "blue", xlim = c(0,300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

summary(mod.exd)
Model fitted: Exponential decay with lower limit at 0 (2 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
d:(Intercept) 0.1167687 0.0021435 54.476 < 2.2e-16 ***
e:(Intercept) 9.3269006 0.2794558 33.375 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error:
0.00182667 (945 degrees of freedom)
1 - pexp(log(1000)*k, d)
[1] 0.005534439
1/d
[1] 7.102853
data %>%
filter(year <= 2010) %>%
pull(prcp) %>%
mean() %>%
log()*k
[1] 23.28785
0.31% change versus 0.442%
Labs:
A Single Years 20 year roll
yr = 1960
data %>%
filter((year >= yr) & (year <= yr + 9)) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
[1] 0.0002201357
plot(x = dens$x, y = fitted(mod), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

yr = 1970
data %>%
filter((year >= yr) & (year <= yr + 9)) %>%
pull(prcp) %>%
hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)
[1] 0.002985969
plot(x = dens$x, y = fitted(mod), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

1 - pexp(log(1000)*k, d)
[1] 0.00457532
mean = 76.84094 sd = 144.0488
Explonential Explainability

Risk Trends
Risks Through out Years
mtr <- data.frame(year = YEAR)
list <- c()
for (i in YEAR){
data %>%
filter((year >= i) & (year <= i + 9)) %>%
pull(prcp) %>%
hist(breaks = 400, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
p <- 1 - pexp(log(1000)*k, d)
list <- c(list, p)
}
NaNs producedNaNs producedNaNs produced
plot(1:length(list), list)

Use smaller breaks
Find the right break numbers
YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22)
lenlist <- c()
kulist <- c()
for(yr in YEAR){
data %>%
filter((year >= yr) & (year <= yr + 19)) %>%
dplyr::select(year, prcp) %>%
pull(prcp) %>%
length() -> n
lenlist <- c(lenlist, n)
data %>%
filter((year >= yr) & (year <= yr + 19)) %>%
dplyr::select(year, prcp) %>%
pull(prcp) %>%
kurtosis() -> ku
kulist <- c(kulist, ku)
}
Square_root = c()
Sturges = c()
Rice = c()
Donna = c()
for(n in lenlist){
sq = ceiling(n^0.5)
sg = ceiling(log(n, 2)) + 1
rc = ceiling(2*n^(1/3))
igh = (6*(n - 2)/(n +1)/(n +3))^0.5
dn = log(abs(1 + kulist[which(lenlist %in% n)]/igh) , 2) + log(n, 2) + 1
Square_root = c(Square_root, sq)
Sturges = c(Sturges, sg)
Rice = c(Rice, rc)
Donna = c(Donna, dn)
}
data.frame( Year = YEAR,
Length = lenlist,
Square_root = Square_root,
Sturges = Sturges,
Rice = Rice,
Donna = Donna
) -> meth
meth
Work Flow
In a 10 years interval compute risks of precipitation higher than 10cm (prcp > 1000) by using different breaks of histogram from usgin 250 to 1000
Choosing bin have many choices actually
YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22, 30, 40, 50, 60, 70, 80, 90)
for(n in BREAKS)
data %>%
filter((year >= 1960) & (year <= 1969)) %>%
pull(prcp) %>%
hist(breaks = n)











kurtosis(data$prcp)
[1] 31.3595
list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
klist <- c()
for(n in BREAKS){
for(i in YEAR){
data %>%
filter((year >= i) & (year <= i + 19)) %>%
pull(prcp) %>%
hist(breaks = n, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$density
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
p <- 1 - pexp(log(1500)*k, d)
list <- c(list, p)
listB <- c(listB, n)
Year <- c(Year, i)
dlist <- c(dlist, d)
elist <- c(elist, e)
klist <- c(klist, k)
}}
mtr <- data.frame(year = Year,
Breaks = listB,
Probability = list,
parameter_d = dlist,
parameter_e = elist,
parameter_k = klist)
Visulisation
ggplotly(g)
n too large, allowed maximum for palette Blues is 9
Returning the palette you asked for with that many colors
mtr %>%
mutate(year = as.factor(year)) %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(y = parameter_d, x = parameter_e, color = year)) +
geom_jitter(alpha = 0.8) +
scale_color_brewer(palette = "Reds") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
) +
xlab("parameter e: skewness of the curve") +
ylab("parameter d: height of the curve")
Parameter d influence the height of pdf, parameter e influence the lenghth of pdf… this suggest in general, precipitation distribution cure are getting flatter and
plot_ly(x= mtr$parameter_e, y= mtr$parameter_d, z=mtr$Breaks, type="scatter3d", mode="markers", color=mtr$year)
mtr %>%
dplyr::select(year, Breaks, parameter_d, parameter_e) %>%
mutate(Breaks = as.factor(Breaks)) %>%
melt(id = c("year", "Breaks"), variable.name = "parameter") %>%
ggplot(aes(x = year, y = value)) +
geom_line(aes(x = year, y = value, linetype = parameter, color = Breaks)) +
scale_color_brewer(palette = "Purples") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
)
mtr %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(x = year, y = parameter_k, color = Breaks)) +
geom_point() +
geom_line() +
scale_color_brewer(palette = "Blues") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
)
calculate different
list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
for(n in BREAKS){
for(i in YEAR){
data %>%
filter((year >= i) & (year <= i + 9)) %>%
pull(prcp) %>%
hist(breaks = n, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids
y <- hist$counts/sum
data.frame(x, y) -> dens
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
p <- 1 - pexp(log(600)*k, d)
list <- c(list, p)
listB <- c(listB, n)
Year <- c(Year, i)
dlist <- c(dlist, d)
elist <- c(elist, e)
}}
mtr <- data.frame(year = Year,
Breaks = listB,
Probability = list,
parameter_d = dlist,
parameter_e = elist)
mtr %>%
mutate(Breaks = as.factor(Breaks)) %>%
ggplot(aes(x = year, y = Probability, color = Breaks)) +
geom_point() +
geom_line() +
scale_color_brewer(palette = "Blues") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
)
g <- df %>%
ggplot(aes(x = date, y = prcp)) +
geom_line(color = "blue") +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#fafafa"),
plot.background = element_rect(fill = "white")
) +
geom_line(aes(x = date, y = rep(1500, length(date))), color = "lightblue") +
ylab("Precipitation (in tenth of mm")
ggplotly(g)
data %>%
filter(prcp >= 1000)
mtr %>%
filter(Breaks == 21) %>%
left_join(meth, by = c("year" = "Year")) %>%
select(year, Probability, Length) %>%
mutate(days = Probability*Length)
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKRGF0YSBmcm9tIDE5NjAgdG8gMjAxMCByYWluaW5nIGRheXMsIHRoZXJlIGFyZSA0MTY4IG9ic2VydmF0aW9ucyAxMyBkYXlzIG9mIHdoaWNoIGFyZSByYWluIGV4Y2VlZHMgMTAwMCAoMTBjbSkKCkNoYW5nZSBvZiBleGNlZWRzIDEwMDAgaXMgMC4wMDMxMTkwMDIgMC4zMSUgY2hhbmdlCgpgYGB7cn0KCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHJlYWRyKQpsaWJyYXJ5KGRyYykKbGlicmFyeShyZXNoYXBlMikKCiNpbnN0YWxsLnBhY2thZ2VzKCJtb21lbnRzIikgICAgICAgICAgICAgICAgCmxpYnJhcnkoIm1vbWVudHMiKSAjIGZvciBza2V3bmVzcyBvZiBkYXRhCmBgYAoKIyBEYXRhCgpgYGB7cn0KcmVhZC5jc3YoIi4vcmF3L3ByY3AuY3N2IikgJT4lCiAgbXV0YXRlKGRhdGUgPSBhcy5EYXRlKGRhdGUpKSAlPiUKICBtdXRhdGUoeWVhciA9IGZvcm1hdChkYXRlLCAiJVkiKSkgJT4lCiAgZmlsdGVyKHByY3AgIT0gMCkgJT4lCiAgZHBseXI6OnNlbGVjdChkYXRlLCB5ZWFyLCBwcmNwKSAtPiBkYXRhCgpkYXRhICU+JQogIGZpbHRlcih5ZWFyIDw9IDIwMjApICU+JQogIHB1bGwocHJjcCkgJT4lCiAgaGlzdChicmVha3MgPSAxMDAwLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKYGBgCgojIE1vZGxpbmcgRGlzdHJpYnV0aW9uCgpVc2UgZGF0YSBmcm9tIDE5NjAgLSAyMDEwIGZvciBtb2RlbGluZyBwcmVjaXBpdGF0aW9uCgpgYGB7cn0KaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQp4IDwtIGhpc3QkbWlkcyAKeSA8LSBoaXN0JGNvdW50cy9zdW0gCgpkYXRhLmZyYW1lKHgsIHkpIC0+IGRlbnMgCmBgYAoKYGBge3J9Cm1vZGVsIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKc3VtbWFyeShtb2RlbCkKZCA9IDAuMTQwNzg4NQplID0gMS4zMjkyNDE0CgprIDwtIDEvZS9kCmBgYAoKYGBge3J9CnBsb3QoeCA9IGRlbnMkeCwgeSA9IGZpdHRlZChtb2RlbCksIHR5cGUgPSAibCIsIGNvbCA9ICJibHVlIiwgeGxpbSA9IGMoMCwgMzAwKSkKbGluZXMoeCA9IGRlbnMkeCwgeSA9IGRlbnMkeSwgdHlwZSA9ICJwIiwgY29sID0gYWxwaGEoImJsYWNrIiwgMC41KSkKYGBgCgpgYGB7cn0KbW9kLmV4ZCA8LSBkcm0oeSB+IHgsIGZjdCA9IEVYRC4yKCksIGRhdGEgPSBkZW5zKQpgYGAKCmBgYHtyfQpwbG90KHggPSBkZW5zJHgsIHkgPSBmaXR0ZWQobW9kLmV4ZCksIHR5cGUgPSAibCIsIGNvbCA9ICJibHVlIiwgeGxpbSA9IGMoMCwzMDApKQpsaW5lcyh4ID0gZGVucyR4LCB5ID0gZGVucyR5LCB0eXBlID0gInAiLCBjb2wgPSBhbHBoYSgiYmxhY2siLCAwLjUpKQpgYGAKCmBgYHtyfQpzdW1tYXJ5KG1vZC5leGQpCmBgYAoKYGBge3J9CjEgLSBwZXhwKGxvZygxMDAwKSprLCBkKQoKMS9kCgpkYXRhICU+JQogIGZpbHRlcih5ZWFyIDw9IDIwMTApICU+JQogIHB1bGwocHJjcCkgJT4lCiAgbWVhbigpICU+JQogIGxvZygpKmsKYGBgCgowLjMxJSBjaGFuZ2UgdmVyc3VzIDAuNDQyJQoKIyBMYWJzOgoKIyMgQSBTaW5nbGUgWWVhcnMgMjAgeWVhciByb2xsCgpgYGB7cn0KeXIgPSAxOTYwCmRhdGEgJT4lCiAgZmlsdGVyKCh5ZWFyID49IHlyKSAmICh5ZWFyIDw9IHlyICsgOSkpICU+JQogIHB1bGwocHJjcCkgJT4lCiAgaGlzdChicmVha3MgPSAxMDAwLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQp4IDwtIGhpc3QkbWlkcyAKeSA8LSBoaXN0JGNvdW50cy9zdW0gCmRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKZCA8LSBtb2QkcGFybU1hdFsxXQplIDwtIG1vZCRwYXJtTWF0WzJdCmsgPC0gMS9lL2QKMSAtIHBleHAobG9nKDEwMDApKmssIGQpCgpwbG90KHggPSBkZW5zJHgsIHkgPSBmaXR0ZWQobW9kKSwgdHlwZSA9ICJsIiwgY29sID0gImJsdWUiLCB4bGltID0gYygwLCAzMDApKQpsaW5lcyh4ID0gZGVucyR4LCB5ID0gZGVucyR5LCB0eXBlID0gInAiLCBjb2wgPSBhbHBoYSgiYmxhY2siLCAwLjUpKQpgYGAKCmBgYHtyfQp5ciA9IDE5NzAKZGF0YSAlPiUKICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyA5KSkgJT4lCiAgcHVsbChwcmNwKSAlPiUKICBoaXN0KGJyZWFrcyA9IDEwMDAsIHBsb3QgPSBGQUxTRSkgLT4gaGlzdApoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCnggPC0gaGlzdCRtaWRzIAp5IDwtIGhpc3QkY291bnRzL3N1bSAKZGF0YS5mcmFtZSh4LCB5KSAtPiBkZW5zIAptb2QgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQpkIDwtIG1vZCRwYXJtTWF0WzFdCmUgPC0gbW9kJHBhcm1NYXRbMl0KayA8LSAxL2UvZAoxIC0gcGV4cChsb2coMTAwMCkqaywgZCkKCnBsb3QoeCA9IGRlbnMkeCwgeSA9IGZpdHRlZChtb2QpLCB0eXBlID0gImwiLCBjb2wgPSAiYmx1ZSIsIHhsaW0gPSBjKDAsIDMwMCkpCmxpbmVzKHggPSBkZW5zJHgsIHkgPSBkZW5zJHksIHR5cGUgPSAicCIsIGNvbCA9IGFscGhhKCJibGFjayIsIDAuNSkpCmBgYAoKYGBge3J9CnlyID0gMTk4MApkYXRhICU+JQogIGZpbHRlcigoeWVhciA+PSB5cikgJiAoeWVhciA8PSB5ciArIDkpKSAlPiUKICBwdWxsKHByY3ApICU+JQogIGhpc3QoYnJlYWtzID0gMTAwMCwgcGxvdCA9IEZBTFNFKSAtPiBoaXN0Cmhpc3QkY291bnRzICU+JSBzdW0oKSAtPiBzdW0KeCA8LSBoaXN0JG1pZHMgCnkgPC0gaGlzdCRkZW5zaXR5IApkYXRhLmZyYW1lKHgsIHkpIC0+IGRlbnMgCm1vZCA8LSBkcm0oZGF0YSA9IGRlbnMsIHkgfiBsb2coeCksIGZjdCA9IEVYRC4yKCkpCmQgPC0gbW9kJHBhcm1NYXRbMV0KZSA8LSBtb2QkcGFybU1hdFsyXQprIDwtIDEvZS9kCjEgLSBwZXhwKGxvZygxMDAwKSprLCBkKQoKZGF0YSAlPiUKICBmaWx0ZXIoKHllYXIgPj0geXIpICYgKHllYXIgPD0geXIgKyA5KSkgJT4lCiAgbXV0YXRlKHRyYW5zID0gbG9nKHByY3ApKmspCmBgYAptZWFuID0gNzYuODQwOTQKc2QgPSAxNDQuMDQ4OAoKIyMgRXhwbG9uZW50aWFsIEV4cGxhaW5hYmlsaXR5CmBgYHtyfQpkYXRhICU+JQogIHB1bGwocHJjcCkgJT4lCiAgaGlzdChicmVha3MgPSAxMDAwLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKaGlzdCRjb3VudHMgJT4lIHN1bSgpIC0+IHN1bQp4IDwtIGhpc3QkbWlkcyAKeSA8LSBoaXN0JGRlbnNpdHkgCmRhdGEuZnJhbWUoeCwgeSkgLT4gZGVucyAKbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKZCA8LSBtb2QkcGFybU1hdFsxXQplIDwtIG1vZCRwYXJtTWF0WzJdCmsgPC0gMS9lL2QKCiMganVzdCB2ZXJpZnkgdGhlIGZ1bmN0aW9uIHdlIHVzZSB3b3JrZWQKZ2dwbG90KGRhdGEuZnJhbWUoeCwgeSksIGFlcyh4LCB5KSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX2Z1bmN0aW9uKGZ1biA9IH4gZCpleHAoLWxvZygueCkvZSksIGNvbG9yID0gInJlZCIpICsgZ2VvbV9saW5lKGRhdGEgPSBkYXRhLmZyYW1lKHgsIGZpdHRlZChtb2QpKSwgYWVzKHgsIGZpdHRlZC5tb2QuKSwgY29sb3IgPSAnYmx1ZScpCgojIGZpdHRlZCBmdW5jdGlvbiAocmVkKSArIGZpdGVkIHZhbHVlIChibHVlKQpnZ3Bsb3QoZGF0YS5mcmFtZSh4LCB5KSwgYWVzKHgsIHkpKSArIGdlb21fcG9pbnQoKSArIGdlb21fZnVuY3Rpb24oZnVuID0gfiBkKmV4cCgtZCpsb2coLngpKmspLCBjb2xvciA9ICJyZWQiKSArIGdlb21fbGluZShkYXRhID0gZGF0YS5mcmFtZSh4LCBmaXR0ZWQobW9kKSksIGFlcyh4LCBmaXR0ZWQubW9kLiksIGNvbG9yID0gJ2JsdWUnKQoKZ2dwbG90KGRhdGEuZnJhbWUoeCwgeSksIGFlcyh4LCB5KSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX2Z1bmN0aW9uKGZ1biA9IH4gZCpleHAoLWQqbG9nKC54KSprKSwgY29sb3IgPSAicmVkIikgKyBnZW9tX2xpbmUoZGF0YSA9IGRhdGEuZnJhbWUoeCwgZml0dGVkKG1vZCkpLCBhZXMoeCwgZml0dGVkLm1vZC4pLCBjb2xvciA9ICdibHVlJykKCmdncGxvdChkYXRhLmZyYW1lKHgsIHkpLCBhZXMoeCwgeSkpICsgZ2VvbV9wb2ludCgpICsgZ2VvbV9mdW5jdGlvbihmdW4gPSB+IGQqZXhwKC1kKmxvZygueCkqayksIGNvbG9yID0gInJlZCIpICsgZ2VvbV9saW5lKGRhdGEgPSBkYXRhLmZyYW1lKHgsIGZpdHRlZChtb2QpKSwgYWVzKHgsIGZpdHRlZC5tb2QuKSwgY29sb3IgPSAnYmx1ZScpCgpgYGAKCgojIyBSaXNrIFRyZW5kcwoKIyMgUmlza3MgVGhyb3VnaCBvdXQgWWVhcnMKCmBgYHtyfQptdHIgPC0gZGF0YS5mcmFtZSh5ZWFyID0gWUVBUikKbGlzdCA8LSBjKCkKCmZvciAoaSBpbiBZRUFSKXsKICAgIGRhdGEgJT4lCiAgICAgIGZpbHRlcigoeWVhciA+PSBpKSAmICh5ZWFyIDw9IGkgKyA5KSkgJT4lCiAgICAgIHB1bGwocHJjcCkgJT4lCiAgICAgIGhpc3QoYnJlYWtzID0gNDAwLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKICAgIGhpc3QkY291bnRzICU+JSBzdW0oKSAtPiBzdW0KICAgIHggPC0gaGlzdCRtaWRzIAogICAgeSA8LSBoaXN0JGNvdW50cy9zdW0gCiAgICBkYXRhLmZyYW1lKHgsIHkpIC0+IGRlbnMgCiAgICBtb2QgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQogICAgZCA8LSBtb2QkcGFybU1hdFsxXQogICAgZSA8LSBtb2QkcGFybU1hdFsyXQogICAgayA8LSAxL2UvZAogICAgcCA8LSAxIC0gcGV4cChsb2coMTAwMCkqaywgZCkKICAgIGxpc3QgPC0gYyhsaXN0LCBwKQp9CnBsb3QoMTpsZW5ndGgobGlzdCksIGxpc3QpCmBgYAoKYGBge3J9CgpgYGAKCiMjIFVzZSBzbWFsbGVyIGJyZWFrcwoKRmluZCB0aGUgcmlnaHQgYnJlYWsgbnVtYmVycwoKYGBge3J9CllFQVIgPC0gYygxOTYwLCAxOTgwLCAyMDAwKQpCUkVBS1MgPC0gYygxOSwgMjAsIDIxLCAyMikKYGBgCgpgYGB7cn0KbGVubGlzdCA8LSBjKCkKa3VsaXN0IDwtIGMoKQpmb3IoeXIgaW4gWUVBUil7CmRhdGEgJT4lCiAgICAgIGZpbHRlcigoeWVhciA+PSB5cikgJiAoeWVhciA8PSB5ciArIDE5KSkgJT4lCiAgICAgIGRwbHlyOjpzZWxlY3QoeWVhciwgcHJjcCkgJT4lIAogICAgICBwdWxsKHByY3ApICU+JQogICAgICBsZW5ndGgoKSAtPiBuCiAgbGVubGlzdCA8LSBjKGxlbmxpc3QsIG4pCmRhdGEgJT4lCiAgICAgIGZpbHRlcigoeWVhciA+PSB5cikgJiAoeWVhciA8PSB5ciArIDE5KSkgJT4lCiAgICAgIGRwbHlyOjpzZWxlY3QoeWVhciwgcHJjcCkgJT4lIAogICAgICBwdWxsKHByY3ApICU+JQogICAgICBrdXJ0b3NpcygpIC0+IGt1CiAga3VsaXN0IDwtIGMoa3VsaXN0LCBrdSkKfQpTcXVhcmVfcm9vdCA9IGMoKQpTdHVyZ2VzID0gYygpClJpY2UgPSBjKCkKRG9ubmEgPSBjKCkKZm9yKG4gaW4gbGVubGlzdCl7CiAgc3EgPSBjZWlsaW5nKG5eMC41KQogIHNnID0gY2VpbGluZyhsb2cobiwgMikpICsgMQogIHJjID0gY2VpbGluZygyKm5eKDEvMykpCiAgaWdoID0gKDYqKG4gLSAyKS8obiArMSkvKG4gKzMpKV4wLjUKICBkbiA9IGxvZyhhYnMoMSArIGt1bGlzdFt3aGljaChsZW5saXN0ICVpbiUgbildL2lnaCkgLCAyKSArIGxvZyhuLCAyKSArIDEKU3F1YXJlX3Jvb3QgPSBjKFNxdWFyZV9yb290LCBzcSkKU3R1cmdlcyA9IGMoU3R1cmdlcywgc2cpClJpY2UgPSBjKFJpY2UsIHJjKQpEb25uYSA9IGMoRG9ubmEsIGRuKQp9CmRhdGEuZnJhbWUoIFllYXIgPSBZRUFSLAogICAgICAgICAgTGVuZ3RoID0gbGVubGlzdCwKICAgICAgICAgICBTcXVhcmVfcm9vdCA9IFNxdWFyZV9yb290LCAKICAgICAgICAgICBTdHVyZ2VzID0gU3R1cmdlcywKICAgICAgICAgICBSaWNlID0gUmljZSwgCiAgICAgICAgICAgRG9ubmEgPSBEb25uYQogICAgICAgICAgICkgLT4gbWV0aApgYGAKCmBgYHtyfQptZXRoCmBgYAoKIyBXb3JrIEZsb3cKCkluIGEgMTAgeWVhcnMgaW50ZXJ2YWwgY29tcHV0ZSByaXNrcyBvZiBwcmVjaXBpdGF0aW9uIGhpZ2hlciB0aGFuIDEwY20gKHByY3AgXD4gMTAwMCkgYnkgdXNpbmcgZGlmZmVyZW50IGJyZWFrcyBvZiBoaXN0b2dyYW0gZnJvbSB1c2dpbiAyNTAgdG8gMTAwMAoKQ2hvb3NpbmcgYmluIGhhdmUgbWFueSBjaG9pY2VzIGFjdHVhbGx5CgpgYGB7cn0KWUVBUiA8LSBjKDE5NjAsIDE5ODAsIDIwMDApCkJSRUFLUyA8LSBjKDE5LCAyMCwgMjEsIDIyLCAzMCwgNDAsIDUwLCA2MCwgNzAsIDgwLCA5MCkKYGBgCgpgYGB7cn0KZm9yKG4gaW4gQlJFQUtTKQogIGRhdGEgJT4lCiAgICAgIGZpbHRlcigoeWVhciA+PSAxOTYwKSAmICh5ZWFyIDw9IDE5NjkpKSAlPiUKICAgICAgcHVsbChwcmNwKSAlPiUKICAgICAgaGlzdChicmVha3MgPSBuKQpgYGAKCmBgYHtyfQprdXJ0b3NpcyhkYXRhJHByY3ApCmBgYAoKYGBge3J9Cmxpc3QgPC0gYygpCmxpc3RCIDwtIGMoKQpZZWFyIDwtIGMoKQpkbGlzdCA8LSBjKCkKZWxpc3QgPC0gYygpCmtsaXN0IDwtIGMoKQpmb3IobiBpbiBCUkVBS1MpewogIGZvcihpIGluIFlFQVIpewogICAgZGF0YSAlPiUKICAgICAgZmlsdGVyKCh5ZWFyID49IGkpICYgKHllYXIgPD0gaSArIDE5KSkgJT4lCiAgICAgIHB1bGwocHJjcCkgJT4lCiAgICAgIGhpc3QoYnJlYWtzID0gbiwgcGxvdCA9IEZBTFNFKSAtPiBoaXN0CiAgICBoaXN0JGNvdW50cyAlPiUgc3VtKCkgLT4gc3VtCiAgICB4IDwtIGhpc3QkbWlkcyAKICAgIHkgPC0gaGlzdCRkZW5zaXR5IAogICAgZGF0YS5mcmFtZSh4LCB5KSAtPiBkZW5zIAogICAgbW9kIDwtIGRybShkYXRhID0gZGVucywgeSB+IGxvZyh4KSwgZmN0ID0gRVhELjIoKSkKICAgIGQgPC0gbW9kJHBhcm1NYXRbMV0KICAgIGUgPC0gbW9kJHBhcm1NYXRbMl0KICAgIGsgPC0gMS9lL2QKICAgIHAgPC0gMSAtIHBleHAobG9nKDE1MDApKmssIGQpCiAgICBsaXN0IDwtIGMobGlzdCwgcCkKICAgIGxpc3RCIDwtIGMobGlzdEIsIG4pCiAgICBZZWFyIDwtIGMoWWVhciwgaSkKICAgIGRsaXN0IDwtIGMoZGxpc3QsIGQpCiAgICBlbGlzdCA8LSBjKGVsaXN0LCBlKQogICAga2xpc3QgPC0gYyhrbGlzdCwgaykKICB9fQptdHIgPC0gZGF0YS5mcmFtZSh5ZWFyID0gWWVhciwKICAgICAgICAgICAgICAgICAgICBCcmVha3MgPSBsaXN0QiwgCiAgICAgICAgICAgICAgICAgICAgICBQcm9iYWJpbGl0eSA9IGxpc3QsIAogICAgICAgICAgICAgICAgICBwYXJhbWV0ZXJfZCA9IGRsaXN0LAogICAgICAgICAgICAgICAgICBwYXJhbWV0ZXJfZSA9IGVsaXN0LAogICAgICAgICAgICAgICAgICBwYXJhbWV0ZXJfayA9IGtsaXN0KQpgYGAKCiMgVmlzdWxpc2F0aW9uCgpgYGB7cn0KbGlicmFyeSh3ZXNhbmRlcnNvbikKbXRyICU+JQogIG11dGF0ZShCcmVha3MgPSBhcy5mYWN0b3IoQnJlYWtzKSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0geWVhciwgeSA9IFByb2JhYmlsaXR5LCBjb2xvciA9IEJyZWFrcykpICsgCiAgICBnZW9tX3BvaW50KCkgKyAKICAgIGdlb21fbGluZSgpICsgCiAgICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZSA9ICJCbHVlcyIpICsKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKSArIAogICAgeWxhYigiUmlza3Mgb2YgUmFpbiBFeGNlZWRzIDEwY20iKSArIAogICAgZ2d0aXRsZSgiUmlza3Mgb2YgUmFpbiBFeGNlZWRpbmcgMTUgY20gaW5jcmVhc2UgZXZleSAyIGRlY2FkZXMiKSArCiAgICB4bGFiKCJUd28gRGVjYWRlcyBUaW1lIEZyb20gXyAiKSAtPiBnCmxpYnJhcnkocGxvdGx5KQpnZ3Bsb3RseShnKQpgYGAKCmBgYHtyfQptdHIgJT4lCiAgbXV0YXRlKHllYXIgPSBhcy5mYWN0b3IoeWVhcikpICU+JQogIG11dGF0ZShCcmVha3MgPSBhcy5mYWN0b3IoQnJlYWtzKSkgJT4lCiAgZ2dwbG90KGFlcyh5ID0gcGFyYW1ldGVyX2QsIHggPSBwYXJhbWV0ZXJfZSwgY29sb3IgPSB5ZWFyKSkgKyAKICAgIGdlb21faml0dGVyKGFscGhhID0gMC44KSArIAogICAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGUgPSAiUmVkcyIpICsKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKSArIAogICAgeGxhYigicGFyYW1ldGVyIGU6IHNrZXduZXNzIG9mIHRoZSBjdXJ2ZSIpICsgCiAgICB5bGFiKCJwYXJhbWV0ZXIgZDogaGVpZ2h0IG9mIHRoZSBjdXJ2ZSIpIAogICAgCmBgYAoKUGFyYW1ldGVyIGQgaW5mbHVlbmNlIHRoZSBoZWlnaHQgb2YgcGRmLCBwYXJhbWV0ZXIgZSBpbmZsdWVuY2UgdGhlIGxlbmdodGggb2YgcGRmLi4uIHRoaXMgc3VnZ2VzdCBpbiBnZW5lcmFsLCBwcmVjaXBpdGF0aW9uIGRpc3RyaWJ1dGlvbiBjdXJlIGFyZSBnZXR0aW5nIGZsYXR0ZXIgYW5kCgpgYGB7cn0KcGxvdF9seSh4PSBtdHIkcGFyYW1ldGVyX2UsIHk9IG10ciRwYXJhbWV0ZXJfZCwgej1tdHIkQnJlYWtzLCB0eXBlPSJzY2F0dGVyM2QiLCBtb2RlPSJtYXJrZXJzIiwgY29sb3I9bXRyJHllYXIpCmBgYAoKYGBge3J9Cm10ciAlPiUKICBkcGx5cjo6c2VsZWN0KHllYXIsIEJyZWFrcywgcGFyYW1ldGVyX2QsIHBhcmFtZXRlcl9lKSAlPiUKICBtdXRhdGUoQnJlYWtzID0gYXMuZmFjdG9yKEJyZWFrcykpICU+JQogIG1lbHQoaWQgPSBjKCJ5ZWFyIiwgIkJyZWFrcyIpLCB2YXJpYWJsZS5uYW1lID0gInBhcmFtZXRlciIpICU+JQogIGdncGxvdChhZXMoeCA9IHllYXIsIHkgPSB2YWx1ZSkpICsgCiAgICBnZW9tX2xpbmUoYWVzKHggPSB5ZWFyLCB5ID0gdmFsdWUsIGxpbmV0eXBlID0gcGFyYW1ldGVyLCBjb2xvciA9IEJyZWFrcykpICsKICAgIHNjYWxlX2NvbG9yX2JyZXdlcihwYWxldHRlID0gIlB1cnBsZXMiKSArCiAgdGhlbWUocGFuZWwuYm9yZGVyID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5tYWpvciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAiI2ZhZmFmYSIpLAogICAgICAgICAgcGxvdC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAid2hpdGUiKQogICAgICAgICkKYGBgCgpgYGB7cn0KbXRyICU+JQogIG11dGF0ZShCcmVha3MgPSBhcy5mYWN0b3IoQnJlYWtzKSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0geWVhciwgeSA9IHBhcmFtZXRlcl9rLCBjb2xvciA9IEJyZWFrcykpICsgCiAgICBnZW9tX3BvaW50KCkgKyAKICAgIGdlb21fbGluZSgpICsgCiAgICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZSA9ICJCbHVlcyIpICsKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKQpgYGAKCmNhbGN1bGF0ZSBkaWZmZXJlbnQKCmBgYHtyfQpsaXN0IDwtIGMoKQpsaXN0QiA8LSBjKCkKWWVhciA8LSBjKCkKZGxpc3QgPC0gYygpCmVsaXN0IDwtIGMoKQpmb3IobiBpbiBCUkVBS1MpewogIGZvcihpIGluIFlFQVIpewogICAgZGF0YSAlPiUKICAgICAgZmlsdGVyKCh5ZWFyID49IGkpICYgKHllYXIgPD0gaSArIDkpKSAlPiUKICAgICAgcHVsbChwcmNwKSAlPiUKICAgICAgaGlzdChicmVha3MgPSBuLCBwbG90ID0gRkFMU0UpIC0+IGhpc3QKICAgIGhpc3QkY291bnRzICU+JSBzdW0oKSAtPiBzdW0KICAgIHggPC0gaGlzdCRtaWRzIAogICAgeSA8LSBoaXN0JGNvdW50cy9zdW0gCiAgICBkYXRhLmZyYW1lKHgsIHkpIC0+IGRlbnMgCiAgICBtb2QgPC0gZHJtKGRhdGEgPSBkZW5zLCB5IH4gbG9nKHgpLCBmY3QgPSBFWEQuMigpKQogICAgZCA8LSBtb2QkcGFybU1hdFsxXQogICAgZSA8LSBtb2QkcGFybU1hdFsyXQogICAgayA8LSAxL2UvZAogICAgcCA8LSAxIC0gcGV4cChsb2coNjAwKSprLCBkKQogICAgbGlzdCA8LSBjKGxpc3QsIHApCiAgICBsaXN0QiA8LSBjKGxpc3RCLCBuKQogICAgWWVhciA8LSBjKFllYXIsIGkpCiAgICBkbGlzdCA8LSBjKGRsaXN0LCBkKQogICAgZWxpc3QgPC0gYyhlbGlzdCwgZSkKICB9fQptdHIgPC0gZGF0YS5mcmFtZSh5ZWFyID0gWWVhciwKICAgICAgICAgICAgICAgICAgICBCcmVha3MgPSBsaXN0QiwgCiAgICAgICAgICAgICAgICAgICAgICBQcm9iYWJpbGl0eSA9IGxpc3QsIAogICAgICAgICAgICAgICAgICBwYXJhbWV0ZXJfZCA9IGRsaXN0LAogICAgICAgICAgICAgICAgICBwYXJhbWV0ZXJfZSA9IGVsaXN0KQoKYGBgCgpgYGB7cn0KbXRyICU+JQogIG11dGF0ZShCcmVha3MgPSBhcy5mYWN0b3IoQnJlYWtzKSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0geWVhciwgeSA9IFByb2JhYmlsaXR5LCBjb2xvciA9IEJyZWFrcykpICsgCiAgICBnZW9tX3BvaW50KCkgKyAKICAgIGdlb21fbGluZSgpICsgCiAgICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZSA9ICJCbHVlcyIpICsKICB0aGVtZShwYW5lbC5ib3JkZXIgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5ncmlkLm1ham9yID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICIjZmFmYWZhIiksCiAgICAgICAgICBwbG90LmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoZmlsbCA9ICJ3aGl0ZSIpCiAgICAgICAgKQpgYGAKCmBgYHtyfQpnIDwtIGRmICU+JQogIGdncGxvdChhZXMoeCA9IGRhdGUsIHkgPSBwcmNwKSkgKyAKICBnZW9tX2xpbmUoY29sb3IgPSAiYmx1ZSIpICsgCiAgdGhlbWUocGFuZWwuYm9yZGVyID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5tYWpvciA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgICAgIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgICBwYW5lbC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAiI2ZhZmFmYSIpLAogICAgICAgICAgcGxvdC5iYWNrZ3JvdW5kID0gZWxlbWVudF9yZWN0KGZpbGwgPSAid2hpdGUiKQogICAgICAgICkgKyAKICBnZW9tX2xpbmUoYWVzKHggPSBkYXRlLCB5ID0gcmVwKDE1MDAsIGxlbmd0aChkYXRlKSkpLCBjb2xvciA9ICJsaWdodGJsdWUiKSArIAogIHlsYWIoIlByZWNpcGl0YXRpb24gKGluIHRlbnRoIG9mIG1tIikKZ2dwbG90bHkoZykKYGBgCgpgYGB7cn0KZGF0YSAlPiUKICBmaWx0ZXIocHJjcCA+PSAxMDAwKQpgYGAKCmBgYHtyfQptdHIgJT4lCiAgZmlsdGVyKEJyZWFrcyA9PSAyMSkgJT4lCiAgbGVmdF9qb2luKG1ldGgsIGJ5ID0gYygieWVhciIgPSAiWWVhciIpKSAlPiUKICBzZWxlY3QoeWVhciwgUHJvYmFiaWxpdHksIExlbmd0aCkgJT4lCiAgbXV0YXRlKGRheXMgPSBQcm9iYWJpbGl0eSpMZW5ndGgpCmBgYAo=